% This matlab script takes the filter parameters as input which is given
% in 'f' and 'choice'(gain for respective frequency components) and order N of the filter gives the
% coefficents of the required filter.

%  This method is implemeted throught the below three stages 
%  1.Define suitable FIR filter
%  2.Read the Input sound and store the data
%  3.Pass through the designed filter, convolute it and obtain the output
%  4.Read the standard signal for the respective genre
%  5.Obtain the delta signal and plot it

% Firstly the sound file is read into an input array using wavread
% function, then it is convoluted with the filter coefficients using the defined function
% conv() 

%  The modified output audio file is played after filtering the input audio
%  file.
%  The Delta signal is plotted.

function [output,fs]=multibandfilter(inputstr,ch,standard)

%###################Define filter #########################################
frequency = [0 0.001875 0.0053125 0.0096875 0.01875 0.03125 0.09375 0.1875 0.375 0.4375 0.5 1];
N=256;

choice=zeros(1,length(ch));
for i=1:length(ch)
choice(i)=(10^(ch(i)/20));
end
%p=choice
filtercoeff = fir2(N,frequency,choice);
[h,w] = freqz(filtercoeff,1,[]);
figure(1);
plot(frequency,choice,w/pi,(abs(h)));
legend('Ideal','fir2 Designed');
title('Comparison of Frequency Response Magnitudes');

% %######read the input sound signal in .wav###################################
%str='file2.wav';
[input,fs]=wavread(inputstr);
%wavplay(input,fs);

input=input(:,1);
figure(2);

subplot(2,1,2);
plot(input);%plotting input signal in time domain
xlabel('time');
ylabel('amplitude');
title('input signal in time domain');  
n=length(input);
[inputf,w]=freqz(input,1,n);
subplot(2,1,1);
plot(w/pi,(abs(inputf)));% plotting input signal in frequency domain 
title('input signal in frequncy domain');
xlabel('normalized frequency');
ylabel('magnitude');
% 
% %###########pass through the filter######################################
%
output=conv(input,filtercoeff);
figure(3);
subplot(2,1,2);
plot(output);
xlabel('time');
ylabel('amplitude');
title('output signal in time domain');
[outputf,w]=freqz(output,1,n);
subplot(2,1,1);
plot(w/pi,(abs(outputf)));
xlabel('normalized frequency');
ylabel('magnitude');
title('output signal in frequency domain');
%hold on;
% 
% %######## play the output#########################
% 
%wavplay(output,fs);

%######## delta signal #########################
if(strcmp(standard,'noname') == 0)
[std,fs]=wavread(standard);
std=std(:,1);
output=output(:,1);
t=length(output)-length(std);
std1=[std ; zeros(t,1)];
y=abs(output)-abs(std1);
figure(4);
plot(y);
title('Delta Signal');
xlabel('Time');
ylabel('Amplitude');
end
